rm(list=ls())

setwd("/Users/mmarbach/Research/ProjectsA/EUSpace/2_EUDim/replication")

library(countrycode)
library(stringr)
library(coda)
library(xtable)
source("fun.R")



#############
# LOAD DATA #
#############

load("./data/for_jags.Rdata")
load("./estimates/out_jags.Rdata")

posterior <- out_jags
cmpemp <- for_jags$cmpemp
phi.select <- for_jags$phi.select

############
# DEVIANCE #
############

deviance <- posterior[,str_detect(colnames(as.matrix(posterior)), paste("deviance",sep="")),] 

plot(deviance)
sapply(deviance, mean)/100

# If there are distinct deviance clusters, select cluster 
# of chains with lowest deviance (=highest posterior probability)
posterior <- mcmc.list(posterior[[5]],posterior[[6]])


###############
# DIAGNOSTICS #
###############

gd1 <- gelman.diag(posterior[,str_detect(colnames(as.matrix(posterior)), paste("K\\[.*,1]",sep="")),])
gd2 <- gelman.diag(posterior[,str_detect(colnames(as.matrix(posterior)), paste("lambda",sep="")),])
gd3 <- gelman.diag(posterior[,str_detect(colnames(as.matrix(posterior)), paste("sigma",sep="")),])

# Drop prior components 
phimcmc1 <- posterior[,str_detect(colnames(as.matrix(posterior)), paste("phi\\[.*,1]",sep="")),]
phimcmc2 <- posterior[,str_detect(colnames(as.matrix(posterior)), paste("phi\\[.*,2]",sep="")),]
K <- ncol(phimcmc1[[1]])
phimcmc1 <- lapply(phimcmc1, function(x) x[,1:(K-12)])
phimcmc2 <- lapply(phimcmc2, function(x) x[,1:(K-12)])
gd4 <- gelman.diag(phimcmc1)
gd5 <- gelman.diag(phimcmc2)

summary(gd1$psrf)
summary(gd2$psrf)
summary(gd3$psrf)
summary(gd4$psrf)
summary(gd5$psrf)


#############
# SUMMARIZE #
#############

# calcuate posterior summaries for "K" (in the paper notation: pi)
mrg.ce <- cmpemp[,c("year","ce", "countryname")]
mrg.ce <- mrg.ce[duplicated(mrg.ce)==FALSE, ]
mrg.ce <- mrg.ce[order(mrg.ce$ce),]

K1 <- posterior[,str_detect(colnames(as.matrix(posterior)), paste("K\\[.*,1]",sep="")),] 
K1mu <- apply(as.matrix(K1), 2, mean) 
K1results <- mrg.ce
mrg.ce$K1mu <- K1mu # mrg.ce must be ordered by 'ce'
K1results$K1mu_base <- K1mu # mrg.ce must be ordered by 'ce'



# Main Figures #
################

cc <- c("Belgium", "France", "Germany", "Italy", "Luxembourg", "Netherlands", 
	"Denmark", "Great Britain", "Ireland", "Greece", "Portugal", "Spain", 
	"Austria", "Finland", "Sweden")
cc_order <- 1:15
dd <- data.frame(countryname=cc, ordering=cc_order)

K1.sampled <- as.matrix(K1)[sample(1:4000,2000),] 

pdf(file=paste("./results/fig_pi_euwest.pdf"), height=7.5, width=4.5, family="Times", pointsize=11)
samp <- subset(mrg.ce, countryname %in% cc) 
samp <- merge(samp, dd, by="countryname") 
samp <- samp[order(samp$ordering),]
fig_countryprob(countries=unique(samp$countryname))
dev.off() 


pdf(file=paste("./results/fig_pi_eueast.pdf"), height=7.5, width=4.5, family="Times", pointsize=11)
samp <- subset(mrg.ce, !( countryname %in% cc) )
fig_countryprob(countries=unique(samp$countryname))
dev.off() 


globaltrend <- trendMaker(post=as.matrix(K1.sampled), 
						  time=mrg.ce$year,
						  newdata=data.frame( time=seq(min(mrg.ce$year), max(mrg.ce$year), by=1))
						  )

pdf(file="./results/fig_pi_avg.pdf",family="Times", width=8.4,height=6.2, pointsize=11)
par(mar=c(2,2,2,0.2), oma=c(0,0,0,0))
plot(1,1,ylim=c(0,1), xlim=c(1945,2010), type="n", ylab="", xlab="",axes=FALSE)
axis(1)
axis(2, at=c(0.2,0.5,0.8))
grid <- seq(0,1,by=0.2)
for(i in 1:length(grid)){
	abline(h=grid[i], lty=1, lwd=0.25, col=grey(0.9))
	}
grid <- seq(1950,2000,by=10)
for(i in 1:length(grid)){
	abline(v=grid[i], lty=1, lwd=0.25, col=grey(0.9))
	}
abline(h=0.5, lty=2)
points(mrg.ce$year, mrg.ce$K1mu, pch=20, col="black")
lines(globaltrend$time, globaltrend$M, col="black",lwd=2)
lines(globaltrend$time, globaltrend$L, col=grey(0.55),lwd=2)
lines(globaltrend$time, globaltrend$U, col=grey(0.55),lwd=2)
dev.off()


	
# Loadings #
############


lambdalabels <- c("Military", "Freedom", "Administration", "Enterprise", "Market", "Protectionism", "Macroeconomics", "Quality of life", "Welfare state", "Traditional morality", "Multiculturalism",  "Labour groups", "Target groups", "Internationalism", "European Union", "National way of life")


lambda1 <- posterior[,str_detect(colnames(as.matrix(posterior)), paste("lambda\\[.*,1]",sep="")),]
lambda2 <- posterior[,str_detect(colnames(as.matrix(posterior)), paste("lambda\\[.*,2]",sep="")),]
lambdamu1 <- as.vector(apply(as.matrix(lambda1), 2, mean))
lambdamu2 <- as.vector(apply(as.matrix(lambda2), 2, mean))
lambdaqu1 <- t(apply(as.matrix(lambda1), 2, quantile, probs=c(0.025, 0.975) ))
lambdaqu2 <- t(apply(as.matrix(lambda2), 2, quantile, probs=c(0.025, 0.975) ))
lambda <- cbind(lambdamu1,lambdamu2)



cat("", file="./results/tab_LambdaResults.tex")
cat( "\\begin{table}[!ht]
\\begin{center}
\\begin{tabular}{lccccc}
\\toprule %%%%%%%%%%%%%%%%%%%%%%%%%
			  & \\multicolumn{2}{c}{LR Dimension} & \\multicolumn{2}{c}{EU Dimension} \\\\
Indicator  & Post. mean $\\lambda$ & BCI & Post. mean $ \\lambda$ & BCI  \\\\
\\midrule %%%%%%%%%%%%%%%%%%%%%%%%% \n"
, append=F, file="./results/tab_LambdaResults.tex")
for(i in 1:16){
	cat(lambdalabels[i], " & " , sprintf("%.2f", lambdamu1[i]), 
											 "& \\scriptsize [" , sprintf("%.2f", lambdaqu1[i,1]) ,", ", sprintf("%.2f", lambdaqu1[i,2]) , "] & ", sprintf("%.2f", lambdamu2[i]), 
											"& \\scriptsize [" , sprintf("%.2f", lambdaqu2[i,1]) ,", ", sprintf("%.2f", lambdaqu2[i,2]) , "] \\\\ \n", 
											sep="", file="./results/tab_LambdaResults.tex", append=T )
	}
cat("\\bottomrule %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\\end{tabular}
\\end{center}
\\end{table}"
, append=T, file="./results/tab_LambdaResults.tex")




# Positions #
#############

K1.sampled <- as.matrix(K1)[sample(1:4000,2000),]
rownames(K1.sampled) <- paste("M_",seq(1,nrow(K1.sampled)),sep="")
K1use <- cbind(K1results,t(K1.sampled))

pitrend <- ddply(K1use, "countryname", function(x){
		countrytrend <- trendMaker(post=t(x[,str_detect(colnames(x),"M_") ]), 
						     time=x$year, 
						      newdata=data.frame( time=x$year)
							)
		return(countrytrend)	
		},.progress = "text")

pitrend$hasEU <- as.numeric(pitrend$M < 0.5)
pitrend <- pitrend[,c("countryname","time", "hasEU")]

phi1 <- as.matrix(posterior[,str_detect(colnames(as.matrix(posterior)), paste("phi\\[.*,1]",sep="")),])
phi1.sampled <- as.matrix(phi1)[sample(1:4000,2000),]
phi2 <- as.matrix(posterior[,str_detect(colnames(as.matrix(posterior)), paste("phi\\[.*,2]",sep="")),])
phi2.sampled <- as.matrix(phi2)[sample(1:4000,2000),]

# merge estimated lr/eu position
phi1 <- as.matrix(posterior[,str_detect(colnames(as.matrix(posterior)), paste("phi\\[.*,1]",sep="")),])
phi2 <- as.matrix(posterior[,str_detect(colnames(as.matrix(posterior)), paste("phi\\[.*,2]",sep="")),])
phi1.mu <- apply(phi1, 2, mean)[1:(ncol(phi1)-8)]
phi2.mu <- apply(phi2, 2, mean)[1:(ncol(phi2)-8)]

phi.mu <- cbind(phi.select, "lr.m1"=phi1.mu, "eu.m1"=phi2.mu)
phi.mu$year <- phi.mu$year * (-1)
toPlot <- phi.mu

toPlot$ord <- NA
toPlot$ord[toPlot$countryname=="Belgium"] <- 1
toPlot$ord[toPlot$countryname=="France"] <- 2
toPlot$ord[toPlot$countryname=="Germany"] <- 3
toPlot$ord[toPlot$countryname=="Italy"] <- 4
toPlot$ord[toPlot$countryname=="Luxembourg"] <- 5
toPlot$ord[toPlot$countryname=="Netherlands"] <- 6
toPlot$ord[toPlot$countryname=="Denmark"] <- 7
toPlot$ord[toPlot$countryname=="Great Britain"] <- 8
toPlot$ord[toPlot$countryname=="Ireland"] <- 9
toPlot$ord[toPlot$countryname=="Greece"] <- 10
toPlot$ord[toPlot$countryname=="Portugal"] <- 11
toPlot$ord[toPlot$countryname=="Spain"] <- 12
toPlot$ord[toPlot$countryname=="Austria"] <- 13
toPlot$ord[toPlot$countryname=="Finland"] <- 14
toPlot$ord[toPlot$countryname=="Sweden"] <- 15


toPlot <- merge(toPlot, pitrend, by.x=c("countryname", "year"), by.y=c("countryname", "time"))

toPlot_ <- toPlot

toPlot <- toPlot[!is.na(toPlot$ord),]

toPlot$sorter <- toPlot$it

filename <- "./results/fig_partyposterior.pdf"
pdf(file=filename, height=7.5, width=4.5,  pointsize=11, family="Times")
par(mfrow=c(5,3),mar=c(2,2,2,0.8), oma=c(0,0,0,0), pch=20)
leftsmallerright <- ddply(toPlot, "ord", function(x){
	plot(0,0,ylim=c(-3,5), xlim=c(-8.5,8), type="n", xlab="", ylab="",  axes=FALSE, main=unique(x$countryname))
  	axis(1); axis(2)
	points(toPlot$lr.m1[toPlot$hasEU==1], toPlot$eu.m1[toPlot$hasEU==1], pch=20, col=gray(0.9))
	x <- x[x$hasEU==1, ]
	if ( nrow(x)==0 ) return(NULL)
  	points(x$lr.m1, x$eu.m1, pch=20, col="black")
  	abline(h=mean(x$eu.m1) - sd(x$eu.m1), lty=3)
  	lr <- phi1.sampled[,x$sorter]
	eu <- phi2.sampled[,x$sorter]
	colnames(lr) <- paste("lr",colnames(lr),sep="_")
	colnames(eu) <- paste("eu",colnames(eu),sep="_")
	draws <- cbind(lr,eu)
  	toPredict <- seq(min(x$lr.m1),max(x$lr.m1),length=30)
  	syst <- apply(draws, 1, function(y){
  		lr_ <- y[str_detect(names(y),"lr_")]
  		eu_ <- y[str_detect(names(y),"eu_")]
  		m <- loess(eu_~lr_, span=10,control = loess.control(surface = "direct"))
  		euhat <- predict(m, newdata=toPredict)
  		return(euhat)
	  	})
  	mu <- apply(syst,1, mean)
  	low <- apply(syst,1, quantile,probs=0.025)
  	hig <- apply(syst,1, quantile,probs=0.975)
  	lines(toPredict, mu,lwd=1.5)
  	lines(toPredict, low,col=grey(0.55),lwd=1.5)
  	lines(toPredict, hig,col=grey(0.55),lwd=1.5)
  	problfrt <- sum(syst[1,] < syst[30,])/nrow(draws)
  	out <- data.frame(countryname=unique(x$countryname),problfrt=problfrt)
  	return(out)
	}, .progress="text")
dev.off()



# Eurosceptic Parties / Taggart #
#################################

cmp <- read.csv("./data/MPDataset_full2010b.csv")
cmp$election <- as.Date(cmp$edate, "%m/%d/%Y")
cmp$year <- format(cmp$election, "%Y")
cmp$share <- (cmp$absseat / cmp$totseats)*100
cmp <- unique(cmp[,c("year", "party", "partyname", "share")])
cmp <- ddply(cmp, c("year","party","partyname"), function(x)  return(share=mean(x$share)) ) 
colnames(cmp)[colnames(cmp)=="V1"] <- "share"

toTab <- toPlot[toPlot$hasEU==1,]
toTab <- merge(cmp, toTab, by.y=c('year', 'party'), by.x=c("year", 'party'), all.y=TRUE, all.x=FALSE)
toTab <- toTab[,c("countryname","year","partyname","party","lr.m1","eu.m1","share")]
eurosceptic <- ddply(toTab, c("countryname", "year"), function(y) {
	y$LR <- ifelse( median(y$lr.m1) > y$lr.m1, "L", "R")
	y <- y[(mean(y$eu.m1) - sd(y$eu.m1)) > y$eu.m1, ]
	return(y)
	})

eurosceptic <- eurosceptic[,c("countryname","year","partyname","party","LR","share")]
print(xtable(as.data.frame(eurosceptic),digits=0), 
	include.rownames=FALSE, 
	file="./results/tab_eurosceptic.tex")

taggart <- read.csv('./data/TaggartID.csv', sep=";")
taggart <- taggart[!(taggart$Country %in% c("Denmark", "Finland")),] 
table( taggart$CMP %in% eurosceptic$party )/length(taggart$CMP)



# Comparisons with Experts #
###########################$

# Reinclude Eastern Europe
toPlot <- toPlot_

ches <- read.csv("./data/1999-2006_ChapelHillSurvey.csv")

parlgov <- read.csv('./data/pargov_party.csv')
convert <- parlgov[,c("cmp", "chess")]
convert <- convert[!is.na(convert$ches), ]
convert <- convert[duplicated(convert)==FALSE,]
colnames(convert) <- c("cmp", "ches")

ches <- merge(convert, ches, by.x="ches", by.y="party_id", all.x=F, all.y=T)
ches$party <- NULL

ches$id <- seq(1,nrow(ches))

ches <- merge(toPlot[toPlot$hasEU==1,], ches, by.x=c("party"), by.y=c("cmp"), suffixes = c("_cmp","_ches") )

ches$diff <- abs(ches$year_cmp - ches$year_ches)
out <- ddply(ches, c("id"), function(x) {
	return(x[which.min(x$diff),])
	} )

cor(out[out$diff < 2,c("lr.m1", "eu.m1", "position", "lrgen")],use="pairwise.complete.obs")

pp <- ddply(out[out$diff < 2, ], "countryname", function(x){
	lr <- cor(x$lr.m1, x$lrgen,use="pairwise.complete.obs")
	eu <- cor(x$eu.m1, x$position,use="pairwise.complete.obs")
	o <- data.frame(eu=eu,lr=lr)
	return(o)
	})

print(xtable(pp,digits=2), 
	include.rownames=FALSE, 
	file="./results/tab_chescorrelations.tex")




